1 Antes de empezar, conozcámonos

Hola, soy Aitor Gonzalez. Yo aprendo tidymodels para que te sea más fácil de aprender.

Soy data scientist o lo que es lo mismo, estadístico que hace más cosas aparte de una regresión logística. Actualmente trabajo en la unidad de facultad de medicina y biología de la Universidad Autónoma de Barcelona.

Este taller contiene mucho texto. Y suele requerir de un par de veces mirárselo para pillarle todo el tranquillo. No pasa nada si no te enteras a la primera. “Dentro de poco” estará subido a Youtube en mi canal Reestimando.

2 ¿Por que tidymodels?

2.1 Un poquito de la historia

La historia comienza con este señor de la foto de la ziquierda. Él es Max Khun. Parece un señor muy majo o muy travieso. Yo personalmente creo más lo segundo viendo su creación.

Comienza por aquí por que este señor fue el desarrollador de la librería caret. caret fue una librería que comenzó en 2008 y que asentó las bases del modelado estadístico en R. Desde su creación, el paquete ha sido citado en más de 7000 publicaciones académicas, teniendo un enorme éxito en la comunidad.

El paquete fue un éxito debido a que antes de él, no había una forma de unificar sintaxis entre los principales paquetes de modelado. Cada modelo estaba en un paquete y todos diferentes entre si, un caos para organizar un flujo de trabajo. caret permitía utilizar con una sintaxis común los modelos de varias librerías. Además también traía varias funciones para hacer automáticamente cosas que son necesarias en todo proceso de modelado, como una función que haga una partición de datos en un conjunto de entrenamiento y uno de prueba (training y test).

ELIPSIS en 2012 aparece la el entorno tidyverse de la mano de Hadley Whickam, la otra súper estrella de R en los últimos años. Tidyverse proporciona un montón de herramientas para poder lidiar con facilidad con la manipulación de datos, convirtiéndose en un package súper popular no solo ne R, sino en toda la ciencia datos. El principal punto es la facilidad de sintaxis para todas sus funciones, haciendo que los usuarios tuvieran una lectura de código muy fácil y fluida, enfocándose 100% en la necesidad y no en la programación.

Encantado con un framework en el que toda la sintaxis está ordenada, en 2019, Max se une al tidyverso, y comienza a crear tidymodels cooperando con el equipo de Whickam. El objetivo era convertir Careta algo más modularizado y que permitiera centrar en el modelado en si mismo, más que infinitas líneas de sintaxis para poder crear modelos de machine learning.

A fecha de este taller, noviembre de 2023, hay otra desarrolladora principal de tidymodels además de Max, Julia Silge.

Esta mujer, no solo es un as de la programación, sino que además se dedica a expandir la filosofía tidy a través de sus redes. Su contenido vale oro y diría que es imperdible si quieres no solo aprender más de tidymodels, sino de las posibilidades del aprendizaje automático en general.

- https://www.youtube.com/@JuliaSilge

- https://juliasilge.com/blog/

Entre algunos de sus trabajos más llamativos están uno que me encanta que es Topic modeling for #TidyTuesday Taylor Swift lyrics es decir, “modelizando los temas que tratan las canciones de Taylor Swift” https://juliasilge.com/blog/taylor-swift/

2.2 Un breve ejemplo

Imaginemos que queremos hacer una regresión lineal. Cogeremos unos datos cualquiera, lo de mtcars.supongamos que queremos hacer una regresión logítica, la sintaxis para eso sería:

Como se pude ver en el ejemplo, la regresión fuera de tidymodels la hacemos con la librería glm. En tidymodels esto se especifica con la función set_engine.

Ahora, veamos que pasaría si queremos hacerla con las liberría h2o

Se va pillando la idea, ¿no?…

tidymodels permite que con cambios mínimos de código puedas cambiar completamente tu flujo de trabajo. Si queremos cambiar un “motor” de modelo, solo tenemos que especificarlo en un solo cambio, y no tener que reescribir todo el código de nuevo. Esto es clave ya que permite jugar mucho a la hora de probar muchos tipos de modelos y ampliar las miras sin aumentar los problemas de código.

Evidentemente, esto no es ni la superficie de tidymodels, aquí solo estamos viendo la isla a lo lejos.

2.3 ¿Que aporta de tidymodels?

  • No más caret: caret está en mantenimiento solo, en el futuro se dará soporte a tidymodels. tidymodels es la apuesta al largo plazo en el entorno de R.
  • Consistencia: Misma sintaxis para todos los procesos.
  • Replicabilidad: Es muy sencillo replicar resultados.
  • Comunicación: Los outputs de tidymodels siguen la lógica de tidyverse.
  • Exportabilidad: Junto a otras librerías, como plumber y vetiver, tidymodels permite fácilmente la exportabilidad de sus modelos a otros entornos en cloud, para poner los modelos a producción.
  • Posibilidades: a fecha de este informe, tidymodels da soporte a 155 modelos, incluyendo modelos para datos censurados y modelos de series temporales.

3 Presentando los datos

Somos unos expertos en aprendizaje automático y nos piden ayuda de un hospital 🏥

Al parecer, nuestra variable, Var_respuesta, es una enfermedad complicada de diagnosticar. Los médicos no son capaces de determinar si un paciente tendrá o no esa enfermedad al cabo de un año. Entonces han hecho un estudio caso-control en el que se miran unas cuantas variables; y al cabo de un año se da un diagnostico definitivo de var_respuesta. En este caso pues, nos interesa una modelización de tipo clasificatoria binaria, es decir de 2 niveles, o tiene la enfermedad o no la tiene. 👍/👎

Además de de nuestra variable Var_respuesta, tenemos otras 25 variables de diversos tipos, de cognición, bioquímicas, genéticas, etc. El dataframe es una simulación de uno real, e incluye valores perdidos y otras características a tener en cuenta en la fase de pre-procesamiento.

Insisto, el propósito es modelizar un clasifiación binaria

Antes de empezar podemos ver un poco su estructura.

#---------------
# Librerias | Carga de datos ----
#---------------
if (!require('pacman') ) {install.packages('pacman')}
pacman::p_load(
  readxl,      # Lectura de los datos
  tidyverse ,  # Acceso al entorno de procesamiento "tidyverse"
  tidymodels,  # Acceso al entorno de modelado "tidymodels"
  agua,        # Modelaje de entorno h2o 
  themis,      # Soporte de la librería "recipes" para tratar desbalanceo
  skimr,       # Descriptiva rápida de datos  
  naniar,      # Resumen y Visualización de datos faltantes
  knitr,       # esta en html para este formato
)

#--- Carga de datos
datos <- read_xlsx('Data/datos.xlsx') %>% 
  filter(!is.na(Var_respuesta) ) %>% 
  mutate(Var_respuesta = as_factor(Var_respuesta))

#--- Primeras 10 filas de los datos.
datos %>% 
  head(.,10) %>% 
  knitr::kable()
Var_respuesta Edad Cognicion_01 Cognicion_02 Cognicion_03 Cognicion_04 Cognicion_05 Bioquimica_01 Bioquimica_02 Bioquimica_03 Bioquimica_04 Bioquimica_05 Bioquimica_06 Bioquimica_07 Bioquimica_08 Bioquimica_09 Bioquimica_10 Bioquimica_11 Escala_01 Escala_02 Escala_03 Genetica_01 Genetica_02 Genetica_03 Genetica_04 Genetica_05 Genetica_06 Genetica_07
0 14 101.82 NA NA NA NA NA NA 5.65 13.65 125.71 NA NA 84.29 1.48 30.35 35.68 0.00000 28.397013 7.363183 NA 0.9665678 NA NA 1.2921746 0.2417377 NA
0 17 80.41 85.53 55.68 124.52 159.70 62.37 NA 5.65 17.83 154.53 NA NA 99.42 1.34 33.99 33.01 138.73119 16.769985 17.764831 0.4550571 NA NA NA NA NA 0.0044324
0 14 91.00 79.76 94.28 140.68 100.07 86.33 26.24 3.18 13.09 135.88 81.53 5.75 109.45 0.95 30.44 34.21 NA 21.382220 7.061147 NA NA NA NA NA NA NA
0 21 86.99 91.33 93.65 115.00 123.01 66.84 10.23 3.48 11.97 132.15 90.66 2.58 68.60 1.04 30.10 34.00 48.83394 7.785552 20.442332 NA NA NA NA NA NA NA
1 33 78.83 104.19 42.68 112.97 179.16 59.51 NA 4.05 13.99 118.33 114.46 3.03 62.87 1.14 22.08 33.50 33.37826 28.165315 28.817018 NA 1.0690024 0.0339304 0.7820054 0.3230982 NA 0.3998781
1 27 93.54 81.16 89.32 115.83 65.07 93.17 32.35 3.50 17.95 132.61 67.61 4.10 76.30 1.11 36.50 32.24 33.45837 7.376139 22.080631 1.3081069 0.5763763 NA 0.5438819 NA 0.4705209 NA
0 16 114.25 NA NA NA NA NA 8.14 5.68 12.25 149.89 NA NA 91.22 1.28 35.86 40.11 NA 18.489371 10.365972 1.5581016 NA NA 1.1597700 NA 0.6242926 0.0579575
0 13 74.80 NA NA NA NA NA 93.32 4.77 8.88 164.33 NA 4.33 84.49 1.07 23.16 27.76 NA 18.749210 9.068513 NA NA NA NA NA NA NA
0 14 66.78 NA NA NA NA NA 54.33 2.56 16.07 161.64 NA NA 99.47 0.92 30.14 34.27 133.14928 18.116769 31.680105 1.6247538 NA 0.9775707 NA 0.9699578 NA 1.1983120
1 17 84.13 100.23 69.62 137.84 91.24 78.76 28.80 4.22 14.38 151.67 NA 5.23 91.19 1.75 18.24 37.76 113.35402 20.327740 21.101902 NA NA 1.1801013 NA NA 1.0568568 NA
#--- Valoración general
skimr::skim(datos)
Data summary
Name datos
Number of rows 285
Number of columns 28
_______________________
Column type frequency:
factor 1
numeric 27
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Var_respuesta 0 1 FALSE 2 0: 227, 1: 58

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Edad 9 0.97 22.33 6.47 8.00 17.00 21.00 27.00 39.00 ▂▇▆▅▂
Cognicion_01 0 1.00 86.51 29.37 0.00 78.50 90.56 103.13 143.84 ▁▁▅▇▂
Cognicion_02 82 0.71 90.22 12.98 55.23 80.74 90.43 98.98 129.15 ▁▆▇▃▁
Cognicion_03 57 0.80 77.62 17.80 35.22 65.91 76.77 88.23 130.09 ▂▇▇▃▁
Cognicion_04 64 0.78 134.85 51.16 15.31 97.17 137.91 172.32 245.38 ▂▅▇▆▂
Cognicion_05 70 0.75 126.23 44.38 53.89 94.13 120.96 155.32 249.33 ▆▇▆▂▁
Bioquimica_01 54 0.81 75.94 14.24 42.80 66.92 74.56 84.50 118.52 ▂▇▇▃▁
Bioquimica_02 43 0.85 36.07 38.35 0.12 11.43 24.47 42.21 250.61 ▇▂▁▁▁
Bioquimica_03 13 0.95 4.05 1.64 1.27 2.90 3.67 4.74 12.68 ▇▇▂▁▁
Bioquimica_04 12 0.96 14.30 1.96 8.88 12.95 14.34 15.54 20.65 ▁▆▇▃▁
Bioquimica_05 18 0.94 140.83 15.09 99.88 129.72 141.38 152.12 187.37 ▁▆▇▃▁
Bioquimica_06 54 0.81 86.32 13.75 55.96 77.21 86.04 95.10 142.74 ▃▇▅▁▁
Bioquimica_07 34 0.88 3.58 1.15 0.83 3.08 3.74 4.30 6.56 ▂▂▇▃▁
Bioquimica_08 16 0.94 84.37 12.35 43.50 76.65 83.96 91.39 126.11 ▁▃▇▃▁
Bioquimica_09 57 0.80 1.13 0.29 0.51 0.91 1.12 1.33 2.09 ▃▇▆▂▁
Bioquimica_10 12 0.96 30.09 3.88 18.24 27.48 30.34 32.51 40.62 ▁▃▇▆▁
Bioquimica_11 14 0.95 33.84 3.66 21.90 31.56 33.71 36.08 43.57 ▁▃▇▆▂
Escala_01 32 0.89 103.18 122.93 0.00 21.23 54.69 137.77 892.17 ▇▂▁▁▁
Escala_02 0 1.00 18.50 8.08 6.17 11.19 18.43 23.58 41.45 ▇▇▇▃▁
Escala_03 0 1.00 18.65 8.24 5.79 12.32 17.76 24.00 46.96 ▇▇▅▂▁
Genetica_01 160 0.44 0.80 0.60 0.00 0.29 0.64 1.23 2.99 ▇▅▃▁▁
Genetica_02 170 0.40 0.78 0.55 0.01 0.33 0.73 1.05 2.27 ▇▇▅▂▂
Genetica_03 178 0.38 0.76 0.56 0.00 0.32 0.67 1.09 2.22 ▇▇▅▂▂
Genetica_04 168 0.41 0.84 0.55 0.01 0.42 0.70 1.20 2.60 ▇▇▅▂▁
Genetica_05 145 0.49 20.19 214.06 0.00 0.84 1.82 3.17 2534.81 ▇▁▁▁▁
Genetica_06 178 0.38 0.84 0.53 0.03 0.37 0.77 1.21 2.36 ▇▇▆▂▁
Genetica_07 137 0.52 0.83 0.59 0.00 0.32 0.72 1.25 2.35 ▇▆▅▃▂
#--- Valores perdidos
naniar::vis_miss(datos)

naniar::gg_miss_upset(datos, nsets = 10, nintersects = 50)

4 ¿De qué se compone tidymodels?

El flujo de trabajo para producir un modelo en tidymodels es el siguiente:

La idea detrás consiste en una sola sintaxis que resuma todo este flujo de trabajo.

  1. Particion de datos
  2. Pre-procesamos los datos de entrenamiento.
  3. Buscamos los híper-parámetros ideales para nuestros datos en el modelo
  4. Ajustamos el modelo.
  5. Comprobamos con las métricas.
  6. Vemos su fnciona con unos datos que el modelo no haya visto.

Como hemos visto en el ejemplo anterior, con cambiar una sola línea línea, podemos cambiar el tipo de modelo. Y con entender unos pocos puntos clave de todos los nexos, podemos llevar a cabo flujos de trabajo muy eficientes con muy poco coste de mantenimiento y poco cambio de sintaxis.

Para la creación de este “cosmos”, tidymodels utiliza las siguientes librerías:

  • rsample: Funciones para crear distintos tipos de remuestreos.
  • recipes: Preprocesador de datos, simplifica el proceso de preparado de datos para cualquier modelo.
  • parsnip: Interfaz para ajustar modelos, hasta 155 tipos a fecha de este html.
  • yarstick: Creación de métricas para la evaluación de un modelo.
  • workflows: Cohesión interna de tidymodels.
  • tune y dials: Creación y gestión de híper-parámetros de cualquier modelo presente en tidymodels.

5 Preparando la partición de datos con Rsampling

Para comenzar a crear un modelo de aprendizaje automático es necesario tener datos, fundamental. Estos datos, suelen ser difíciles de conseguir, requieren tiempo, esfuerzo y dinero. Su manejo debe ser optimizado con el fin de abaratar los costes de recopilarles y mejorar los resultados. Además tener en cuenta que el propósito de un modelo de aprendizaje automático en general es que se pueda usar para poder predecir que va a pasar con unos datos que no tenga, o que tenga en un futuro.

Aquí entre el punto del remuestreo. Crearemos un partición de la muestra en unos datos para “entrenar” el modelo y la otra para “evaluarlo”, para ver si tal y como ha quedado calibrado el modelo es útil o no.

Por supuesto, la calidad de esta partición estará ligada a la calidad de los datos. Si los datos recopilados no son suficientemente representativos, el modelo no predecirá bien cuando traigamos casos que no haya visto, aún incluso si la evaluación del modelo ha sido buena.

Básicamente queremos esto:

La librería dentro de tidymodels que nos dará los conjuntos de datos para hacer esto es rsample.

Algunas de las funciones importantes dentro de este conjunto:

  • initial_split(): Principal función dentro de la librería.Crea un objeto “split” sencillo para hacer una partición de datos
# 🔴 La función principal de Rsample, hace una partición de un conjunto de  🔴
initial_split(datos, prop = 0.6 )
## <Training/Testing/Total>
## <171/114/285>
  • bootstraps(): crea un tibble con n cantidad de objetos “split” y su correspondiente identificador. Estas muestras son de tipo bootstrap, es decir, muestras en las que un individuos puede aparecer más de una vez en el conjunto de datos, tanto en entrenamiento como en evaluación.
rsample::bootstraps(datos,5)
## # Bootstrap sampling 
## # A tibble: 5 × 2
##   splits            id        
##   <list>            <chr>     
## 1 <split [285/99]>  Bootstrap1
## 2 <split [285/103]> Bootstrap2
## 3 <split [285/113]> Bootstrap3
## 4 <split [285/104]> Bootstrap4
## 5 <split [285/113]> Bootstrap5
  • rolling_origin(): Crea un tibble con “n” cantidad de objetos “split” y su correspondiente identificador. En este caso, los remuestreos no son aleatorios y cada partición contienen datos que consecutivos. Es decir, se comienza por la fila 1 y se hacen bloqueen a cantidades iguales. Se usa sobretodo en la partición de datos de tiempo o cuando se está estudiando la creación de una partición optimizada entre el conjunto de entrenamiento y de evaluación.
rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = FALSE)
## # Rolling origin forecast resampling 
## # A tibble: 6 × 2
##   splits          id    
##   <list>          <chr> 
## 1 <split [20/20]> Slice1
## 2 <split [20/20]> Slice2
## 3 <split [20/20]> Slice3
## 4 <split [20/20]> Slice4
## 5 <split [20/20]> Slice5
## 6 <split [20/20]> Slice6
ggpubr::ggarrange(
  rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = FALSE) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile() + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  rolling_origin(datos, initial = 50, assess = 20, skip = 70, cumulative = FALSE) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile()  + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = T) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile()  + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  rolling_origin(datos, initial = 50, assess = 50, skip = 70, cumulative = T) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile()  + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  common.legend = T, legend = 'bottom'
)

5.1 Estratificar

Una de las cosas que pasan habitual a la hora de crear una partición de datos es que, al ser aleatorias, uno de los conjuntos contenga más datos de un tiempo de otro tipo.

Idealmente querríamos que el conjunto de entrenamiento y que el conjunto de evaluación contuvieran la misma proporción de las categorías.

Para ello, en varias de las funciones de rsample se cuenta con el argumento strata que permite hacer esto de forma estratificada para una variable o variables de interés, es este caso, para nuestros datos y nuestra variable conjunta:

set.seed(4147)
# Splits ----

## Sin estratificar 
Split_datos_no_strat <- initial_split(datos, prop = 0.70 )
# Put 3/4 of the data into the training set 
# 🔴 Initial_split con datos estratificados 🔴
## Estratificados 
Split_datos_strat <- initial_split(datos, prop = 0.70, strata = Var_respuesta)
ggpubr::ggarrange(
  ggpubr::ggarrange(
    training(Split_datos_no_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      ungroup() %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
    scale_fill_manual( values = c("azure4", "azure3"))
    ,
    testing(Split_datos_no_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ungroup() %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("azure4", "azure3")),
    common.legend = T, legend = 'bottom'
  ),
  ggpubr::ggarrange(
    training(Split_datos_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      ungroup() %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("#E86EB7", "#108064"))
    ,
    testing(Split_datos_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ungroup() %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("#E86EB7", "#108064")), 
    common.legend = T,
    legend = 'bottom'
  ),
  labels = c('No estratificados', 'Estratificados'), label.x = 0.1
) %>%  
  plot()

Aclaración: los colores elegidos para el gráfico son de mi escena preferida de Spiderman, into the spiderverse. dejad siempre “huevos de Pascua” en los trabajos, hace que las cosas sean más emocionantes :P

5.2 Deshacer la partición

Para esto usaremos las funciones training() y test()

# 🔴 Deshacemos las particiones con las funciones training y testing 🔴

Training_datos <- training(Split_datos_strat)
Testing_datos  <- testing(Split_datos_strat)

knitr::kable(head(Training_datos,10))
Var_respuesta Edad Cognicion_01 Cognicion_02 Cognicion_03 Cognicion_04 Cognicion_05 Bioquimica_01 Bioquimica_02 Bioquimica_03 Bioquimica_04 Bioquimica_05 Bioquimica_06 Bioquimica_07 Bioquimica_08 Bioquimica_09 Bioquimica_10 Bioquimica_11 Escala_01 Escala_02 Escala_03 Genetica_01 Genetica_02 Genetica_03 Genetica_04 Genetica_05 Genetica_06 Genetica_07
0 14 66.78 NA NA NA NA NA 54.33 2.56 16.07 161.64 NA NA 99.47 0.92 30.14 34.27 133.149284 18.116769 31.680105 1.6247538 NA 0.9775707 NA 0.9699578 NA 1.1983120
0 23 108.18 84.49 82.39 183.52 166.31 101.24 3.66 3.18 20.29 148.58 100.23 3.53 81.72 1.11 35.75 34.03 68.943373 12.183308 15.670177 NA 1.797438 NA NA 0.1373531 1.1481697 NA
0 26 0.00 NA NA NA NA NA 129.41 5.99 14.98 145.04 94.53 4.20 78.76 0.87 32.57 26.89 2.898493 11.236520 6.628174 NA NA NA 1.6622978 2.6878823 1.6836672 0.1001570
0 14 52.88 NA NA NA NA NA NA 2.06 13.98 132.94 105.75 NA 82.49 NA 34.73 29.83 99.600374 25.026820 12.949102 0.1247104 NA NA 0.6977428 1.2486785 0.4195379 0.3760827
0 27 109.25 88.04 103.31 143.93 144.97 100.08 8.09 4.49 19.09 187.37 84.47 3.87 91.97 1.26 30.74 34.95 0.000000 16.798451 7.607729 0.8493442 NA 0.8781136 1.9962258 NA NA NA
0 30 105.60 104.09 86.45 149.99 105.02 93.65 7.20 4.51 15.65 145.00 90.13 3.57 90.67 1.05 28.98 31.50 51.087736 6.731747 10.894387 NA 1.349131 0.3128530 NA 1.9916249 1.2268437 NA
0 12 91.73 NA NA NA NA NA 30.29 3.15 13.94 134.92 103.40 3.97 96.92 1.27 21.42 29.26 160.752695 17.851526 18.519326 0.8739840 NA 0.5427362 NA 1.2463796 0.8887141 0.7917460
0 17 91.97 95.35 80.65 172.32 122.78 68.50 40.55 2.77 12.49 145.84 74.00 4.01 86.12 1.16 27.01 33.90 15.319511 14.058768 17.484487 0.2352628 NA 0.3039570 NA NA NA 0.9261596
0 14 69.65 NA NA NA NA NA 55.17 4.87 11.98 165.35 91.39 5.73 52.67 1.21 29.25 32.46 76.930164 22.999587 14.681546 NA NA NA NA NA NA NA
0 34 111.50 58.99 101.17 181.40 119.87 106.90 9.96 3.19 14.47 159.95 102.97 3.98 85.87 0.91 30.43 29.18 7.373436 8.195322 6.718927 1.2862183 NA NA NA NA NA 0.2123626
knitr::kable(head(Testing_datos,10))
Var_respuesta Edad Cognicion_01 Cognicion_02 Cognicion_03 Cognicion_04 Cognicion_05 Bioquimica_01 Bioquimica_02 Bioquimica_03 Bioquimica_04 Bioquimica_05 Bioquimica_06 Bioquimica_07 Bioquimica_08 Bioquimica_09 Bioquimica_10 Bioquimica_11 Escala_01 Escala_02 Escala_03 Genetica_01 Genetica_02 Genetica_03 Genetica_04 Genetica_05 Genetica_06 Genetica_07
0 14 101.82 NA NA NA NA NA NA 5.65 13.65 125.71 NA NA 84.29 1.48 30.35 35.68 0.00000 28.397013 7.363183 NA 0.9665678 NA NA 1.292175 0.2417377 NA
0 17 80.41 85.53 55.68 124.52 159.70 62.37 NA 5.65 17.83 154.53 NA NA 99.42 1.34 33.99 33.01 138.73119 16.769985 17.764831 0.4550571 NA NA NA NA NA 0.0044324
0 14 91.00 79.76 94.28 140.68 100.07 86.33 26.24 3.18 13.09 135.88 81.53 5.75 109.45 0.95 30.44 34.21 NA 21.382220 7.061147 NA NA NA NA NA NA NA
0 21 86.99 91.33 93.65 115.00 123.01 66.84 10.23 3.48 11.97 132.15 90.66 2.58 68.60 1.04 30.10 34.00 48.83394 7.785552 20.442332 NA NA NA NA NA NA NA
0 16 114.25 NA NA NA NA NA 8.14 5.68 12.25 149.89 NA NA 91.22 1.28 35.86 40.11 NA 18.489371 10.365972 1.5581016 NA NA 1.159770 NA 0.6242926 0.0579575
0 13 74.80 NA NA NA NA NA 93.32 4.77 8.88 164.33 NA 4.33 84.49 1.07 23.16 27.76 NA 18.749210 9.068513 NA NA NA NA NA NA NA
1 17 84.13 100.23 69.62 137.84 91.24 78.76 28.80 4.22 14.38 151.67 NA 5.23 91.19 1.75 18.24 37.76 113.35402 20.327740 21.101902 NA NA 1.180101 NA NA 1.0568568 NA
0 31 128.68 NA 99.23 173.48 144.55 101.62 21.34 5.61 14.16 147.96 104.95 4.06 91.39 1.23 31.89 32.58 118.10758 17.200513 14.212995 NA NA NA 0.972422 1.642294 NA 0.5527762
0 20 110.69 88.69 72.28 143.15 193.88 86.83 25.79 4.50 13.12 134.66 68.58 3.56 81.48 1.11 28.27 31.64 225.00452 17.789336 20.238253 NA NA NA NA NA NA NA
0 16 83.73 NA NA NA NA NA 23.32 4.16 12.12 137.10 83.00 5.53 71.56 0.99 32.46 35.09 133.40626 20.677244 15.434748 NA NA NA NA NA NA NA

5.3 folds, creando conjuntos de validación

# 🔴 Creamos Folds 🔴
Folds_training <- vfold_cv(Training_datos, v = 5 ,strata = Var_respuesta)

6 Pre-procesamiento de datos con parsnip

El pre-procesamiento es clave en el proceso de creación de un modelos de machine learning. Cuando queremos hacer un dataframe para modelizar en aprendizaje automático, no viene todo arreglado. Incluso con una recopilación de datos buena puede que necesitemos algunos arreglos. Algunas de las causas de esto pueden ser:

  • Hay modelos que no aceptan valores perdidos. Es necesario que imputemos (rellenar con valores perdidos).
  • Dependiendo de cómo queramos escalar el modelo o a qué pregunta estemos contestando, los datos deberán estar normalizados o no.
  • Re-equilibrar grupos disparejos aumenta la eficiencia de predicción, aunque hay que tener cuidado de cuando usarlo.
  • Reconvertir/transformar variables para maximizar el rendimiento del modelo.

Hacer todo estos procesos puede ser engorroso, pero podemos salir del paso fácilmente gracias por ejemplo de dplyr. Sin embargo ¿Cómo hacemos cuando queramos poner un modelo a grane scala y entren datos a cada momento? ¿Debemos ejecutar el pre-procesamiento una y otra vez?.

La clave para resolver esto está en el paquete recipes. Esta librería trae consigo todo un conjunto de funciones que permiten un pre-procesamiento fino, fino.

6.1 ¿Cómo cocinamos una receta?

Para poder elaborar un flujo completo de pre-procesamiento debemos seguir los siguuientes pasos:

  1. recipe(): la base. Requiere de qué fórmula queremos para el modelo y de qué conjunto de datos partimos. Los datos aquí incorporados serán los datos de entrenamiento
  2. step(): el proceso que queramos aplicar a nuestros datos, cualquier transformación de datos va aquí.
  3. prep(): la receta se prepara, es decir se aplica al conjunto de entrenamiento y se abstrae para un modelado cualquiera
  4. bake(): pone el pre-procesamiento en producción para un conjunto de datos cualquiera.

Hagamos un ejemplo básico: queremos una receta que sencillamente impute los datos faltantes por la media correspondiente a cada variable. Para eso vamos con ejemplo sencillo. La siguiente receta dejará los datos preparados para un modelo de aprendizaje automático, pero sin aplicarle ninguna transformación.

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. step
  # 🟡 En este caso no aplicamos ninguina trasnformación 
  # Step_log
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>% 
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL) %>% 
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()
Cognicion_01 Cognicion_02 Escala_01 Escala_02 Var_respuesta
66.78 NA 133.149284 18.116769 0
108.18 84.49 68.943373 12.183308 0
0.00 NA 2.898493 11.236520 0
52.88 NA 99.600374 25.026820 0
109.25 88.04 0.000000 16.798451 0
105.60 104.09 51.087736 6.731747 0
91.73 NA 160.752695 17.851526 0
91.97 95.35 15.319511 14.058768 0
69.65 NA 76.930164 22.999587 0
111.50 58.99 7.373436 8.195322 0
Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. Step (haz la transofrmación que requieras) 
  # 🟡 Este step será el logaritmo de todas las variables predictoras númericas.
  step_log(all_numeric(), -all_outcomes()) %>% 
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>% 
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL) %>% 
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()
Cognicion_01 Cognicion_02 Escala_01 Escala_02 Var_respuesta
4.201404 NA 4.891471 2.896838 0
4.683796 4.436633 4.233286 2.500067 0
-Inf NA 1.064191 2.419169 0
3.968025 NA 4.601166 3.219948 0
4.693639 4.477791 -Inf 2.821287 0
4.659658 4.645256 3.933544 1.906835 0
4.518849 NA 5.079867 2.882089 0
4.521462 4.557554 2.729127 2.643246 0
4.243483 NA 4.342898 3.135476 0
4.714025 4.077368 1.997884 2.103563 0

¡WALA! Datos pre-procesados y listos.

A partir de aquí en cielo es el límite ya que hay steps para todo, recomiendo leer la librería recipes y ver que tienen disponible. Al final de esta sección dejaremos lista una receta como ejemplo

6.2 Un par de ejemplos más

6.2.1 Imputación

Hacer una imputación simple suele no ser la mejor vía de hacerlo. lo ideal siempre es imputación múltiple. Es decir, una imputación que rellena los valores que faltan con números plausibles derivados de las distribuciones y relaciones entre las variables observadas en el conjunto de datos.

En tidymodels se puede hacer imputación múltiple, y valga la redundancia, de múltiples formas. Por ejemplo, podemos imputar una variable usando una regresión lineal a partir de otras. Hagamos una prueba con la variable Cognicion_02 a partir de Cognicion_01 y Escala_02

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. Step (haz la transofrmación que requieras) 
  # 🟡 Este step será la imputación por regresiópn lineal con otras covariables
  step_impute_linear( Cognicion_02, impute_with = imp_vars(Cognicion_01, Escala_02) ) %>%
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>%
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL) %>%   
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()
Cognicion_01 Cognicion_02 Escala_01 Escala_02 Var_respuesta
66.78 95.30378 133.149284 18.116769 0
108.18 84.49000 68.943373 12.183308 0
0.00 105.71896 2.898493 11.236520 0
52.88 97.79188 99.600374 25.026820 0
109.25 88.04000 0.000000 16.798451 0
105.60 104.09000 51.087736 6.731747 0
91.73 91.30365 160.752695 17.851526 0
91.97 95.35000 15.319511 14.058768 0
69.65 95.03224 76.930164 22.999587 0
111.50 58.99000 7.373436 8.195322 0

Esto va genial cuando conocemos la estructura de los datos y sabemos con certeza la dependencia de unas variables con otras para imputarlas fácilmente, en lugar de confiar ciegamente en algoritmos de imputación.

6.2.2 PCA

Como parte del pre-procesamiento, a veces puede ser necesario convertir ciertas variables a PCA para usarlas como componentes. No entraré en su uso, solo especificar la forma de hacerlo.

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 + Cognicion_02 + Cognicion_03 + Escala_01 + Escala_02  + Escala_03"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. Step (haz la transofrmación que requieras) 
  # 🟡 Este step será la imputación por regresiópn lineal mcon otras covariables
  step_impute_bag(all_numeric()) %>% 
  step_pca( Cognicion_01 , Cognicion_02 , Cognicion_03, num_comp = 2,  id = "pca") %>%
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>%
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  tidy(id = "pca") %>%
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()
terms value component id
Cognicion_01 -0.5964794 PC1 pca
Cognicion_02 -0.6133048 PC1 pca
Cognicion_03 -0.5177543 PC1 pca
Cognicion_01 0.7421220 PC2 pca
Cognicion_02 -0.6671338 PC2 pca
Cognicion_03 -0.0647105 PC2 pca
Cognicion_01 0.3057241 PC3 pca
Cognicion_02 0.4228354 PC3 pca
Cognicion_03 -0.8530785 PC3 pca

6.2.3 Re-remuestrear

ADASYN, o Adaptive Synthetic Sampling, es un algoritmo bastante utilizado cuando tenemos conjuntos de datos desequilibrados.

Los conjuntos de datos desequilibrados se producen cuando la distribución de clases en los datos es desigual, lo que dificulta el entrenamiento de modelos que pueden estar sesgados hacia la clase mayoritaria.

https://www.researchgate.net/publication/224330873_ADASYN_Adaptive_Synthetic_Sampling_Approach_for_Imbalanced_Learning

# 🔴 Libreria que permite instalar algoritmos e sobremuestreo en tidymodels.
pacman::p_load(themis)

Training_datos_adasyn <- Training_datos %>% 
  # 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  #  2. Step (haz la transofrmación que requieras) 
  #  step de  imputación por arboles manteniendo la estrucura de todas las covariables
  step_impute_bag(all_numeric()) %>% 
  # 🔴 step para sobremuestrear la el nivel menor en la variable respuesta.
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  #  3. Prep (asienta la receta y generalizala) 
  prep() %>%
  #  4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL)
# mostramos resultados
Training_datos_adasyn %>% 
  head(.,10) %>% 
  knitr::kable()
Cognicion_01 Cognicion_02 Escala_01 Escala_02 Var_respuesta
66.78 102.05505 133.149284 18.116769 0
108.18 84.49000 68.943373 12.183308 0
0.00 99.58442 2.898493 11.236520 0
52.88 97.03067 99.600374 25.026820 0
109.25 88.04000 0.000000 16.798451 0
105.60 104.09000 51.087736 6.731747 0
91.73 89.86937 160.752695 17.851526 0
91.97 95.35000 15.319511 14.058768 0
69.65 100.57019 76.930164 22.999587 0
111.50 58.99000 7.373436 8.195322 0
ggpubr::ggarrange(
  Training_datos %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      ungroup() %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(
        aes(label = scales::percent(prop_var_respuesta)), 
        position = position_stack(vjust = 0.5), size = 4) +
    scale_fill_manual( values = c("azure4", "azure3"))
    ,
    Training_datos_adasyn %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ungroup() %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(
        aes(label = scales::percent(prop_var_respuesta)),
        position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("azure4", "azure3")),
    common.legend = T, legend = 'bottom'
  ,
  labels = c('Normal', 'Sobre Muestreo\n        Adasyn'), label.x = 0.1
) %>%  
  plot()

6.3 Caso práctico: receta

Por el momento dejaré una receta lista para continuar con el taller.

# Creamos la Fórmula para la receta a partir 
# de nombres de las variables y un poco de magia con paste
Receta_modelo_1_formula <- Training_datos %>% 
  dplyr::select(
   dplyr::matches('Cog'), dplyr::matches('escala')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_1 <- recipe(  
    formula = Receta_modelo_1_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 20% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.1) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_bag(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()
Receta_modelo_1 %>% 
  bake(., new_data = NULL) %>% 
  head(., 10) %>% 
  knitr::kable()
Cognicion_01 Escala_02 Escala_03 Var_respuesta
-0.6426710 -0.0464912 1.5341196 0
0.7604689 -0.7638606 -0.3800329 0
-2.9059967 -0.8783295 -1.4610980 0
-1.1137735 0.7889501 -0.7053657 0
0.7967336 -0.2058789 -1.3439820 0
0.6730268 -1.4229671 -0.9510280 0
0.2029411 -0.0785598 -0.0393877 0
0.2110752 -0.5371131 -0.1631135 0
-0.5454002 0.5438529 -0.4982340 0
0.8729912 -1.2460174 -1.4502475 0

7 ¿Dónde está el modelo? Creando modelos con parnsip

Vale, dejamos atrás la parte del pre-procesamiento y empezamos la parte más machine leaninrg. Modelos. Modelos a la carrera.

La función del paquete parsnip dentro de tidymodels es propoconar una interfaz universal para diferentes métodos de modelado en R.

El ejemplo del principio del taller está basado en esta sección, parnsip es el encargado de que solo con cambiar unas pocas líneas puedas cambiar todo un modelo dentro de un flujo de aprendizaje automático.

A fecha de este taller, tidymodels cuenta con más de 155 modelos, incluyendo supervivencia y series de tiempo. Se pueden encontrar en: https://www.tidymodels.org/find/parsnip/

Supongamos que queremos empezar modelar, queremos usar una regresión logística, un clásico en este campo.Pero dentro de R hay varias librarías que usan la regresión logística, cada una con sus variaciones y reglas. Veamos que tenemos para esto:

# 🔴 Ver los posibles motores de un tipo de modelo
show_engines('logistic_reg') 
## # A tibble: 8 × 2
##   engine    mode          
##   <chr>     <chr>         
## 1 glm       classification
## 2 glmnet    classification
## 3 LiblineaR classification
## 4 spark     classification
## 5 keras     classification
## 6 stan      classification
## 7 brulee    classification
## 8 h2o       classification
logistic_reg(
  engine = "glmnet",
  #Hyper-parametros           
  penalty = NULL, mixture = NULL, mode = 'classification')
## Logistic Regression Model Specification (classification)
## 
## Computational engine: glmnet
logistic_reg(
  engine = "glmnet",
  #Hyper-parametros           
  penalty = 0, mixture = 0, mode = 'classification')
## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = 0
##   mixture = 0
## 
## Computational engine: glmnet
logistic_reg(
  engine = "glmnet",
  #Hyper-parametros           
  penalty = tune(), mixture = tune(), mode = 'classification')
## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = tune()
##   mixture = tune()
## 
## Computational engine: glmnet

7.1 Caso práctico, modelo de Random Forest

Model_RandomForest <- 
  # 🔴 Especificamos el modelo que queremos en este caso un random forest
  rand_forest(
  # 🟡 Ajustamos los hiper-parametros 
    mtry = tune(),  trees = tune(),  min_n = tune()) %>% 
  # 🔴 Ponemos el motor, es decir la librería por la queremos que se ejecute el modelo
  set_engine("ranger", importance = "impurity") %>% 
  # 🔴 Ajustamos el modo del modelo, es decri si queremos regresion o clasifiacion 
  set_mode("classification")

8 Unificando todo: Workflow

El paquete de workflows permite combinar una receta y una especificación de modelo en un único objeto. Esto hace que todo el proceso de modelado, desde el pre-procesamiento de datos hasta el ajuste del modelo, esté en un mismo punto

WF_Receta_modelo_1_Random_forest <- 
  # 🔴 Activamos el workflow 
  workflow() %>% 
  # 🔴 Añadimos la receta
  add_recipe(Receta_modelo_1) %>% 
  # 🔴 Añadimos el modelo
  add_model(Model_RandomForest)

9 Eligiendo métricas con yardstick

El paquete yardstick, en concreto, se utiliza para las métricas de evaluación de modelos. Además

# 🔴 Seleccionamos las métricas que queremos para nuestro modelo
Modelo_Metricas <- metric_set(accuracy, j_index, precision, sensitivity, specificity, roc_auc, f_meas,recall, mcc)

También decir que es es posible crear tu propia métrica. Pero aún no he indagado en esta metodología. Se puede encontrar más en: https://yardstick.tidymodels.org/articles/custom-metrics.html

10 Eligiendo hiperparámetros: tune y dials

Los hiperparámetros son una pieza crucial en los modelos de aprendizaje automático. Estos son parámetros de configuración externos que no pueden aprenderse a partir de los datos, pero que influyen significativamente en el rendimiento del modelo. Un ejemplo sería el número de árboles de un random forest por ejemplo. A diferencia de los parámetros del modelo, que se aprenden durante el entrenamiento, los híperparámetros se establecen antes de que comience el entrenamiento de este.

La selección de los hiperparámetros adecuados es esencial para conseguir un modelo de aprendizaje automático generalizable y con buen rendimiento y es un proceso delicado. Los híper-parámetros desempeñan un papel crucial a la hora de controlar el equilibrio entre la sobreajuste y el infraajuste (overfitting y underfitting). El sobreajuste se produce cuando un modelo funciona bien con los datos de entrenamiento pero no consigue generalizarse a datos nuevos que no se han visto. Por otro lado, el infraajuste se produce cuando el modelo es demasiado simple y no puede captar los patrones subyacentes en los datos. Ajustar los hiperparámetros puede ayudar a encontrar el nivel adecuado de complejidad del modelo, haciendo un buen equilibrio entre la varianza del conjunto de entrenamiento y el de evaluación.

Los paquetes tune y dials son los que nos permiten hacer esta selección y control. Aquí se crea todo un subflujo de trabajo para el ajuste de híper-parámetros.

  1. Creamos unas réplicas del conjunto de entrenamiento, para tener subconjunto de entrenamiento y de validación (que ya hicimos en el apartado de rsample con el objeto ““)
  2. Elegimos un modelo con parsnip y dejamos sus híperparametros en abierto con tune().
  3. Creamos una cuadrícula (grid) con diferentes posibles combinaciones de los híper-parámetros que dejamos en abierto en el paso 2.
  4. Ponemos todo junto con la función tune_grid()

10.1 Montando un espacio de híper-parametros: familia grid_

  • grid_regular(): hará todas las combinaciones entre los rangos de la variable.
  • grid_random(): probará hiper parámetros aleatoriamente.
  • grid_max_entropy(): propondrá una combinacion de hiperpárametros que garanticen que se cubrirá todo el espectro, con fin de evitar mínimos locales en la conversión de algoritmos.
  • grid_latin_hypercube(): propondrá una matriz.

grid_regular() hará las combinaciones acorde al rango que le hayamos dado. Por defecto hace la combinación entre el mínimo, el medio y el máximo de cada valor.

grid_regular(
  mtry(range  = c(2, 10)),
  min_n(range = c(2, 10)),
  trees(range = c(2, 10)))
## # A tibble: 27 × 3
##     mtry min_n trees
##    <int> <int> <int>
##  1     2     2     2
##  2     6     2     2
##  3    10     2     2
##  4     2     6     2
##  5     6     6     2
##  6    10     6     2
##  7     2    10     2
##  8     6    10     2
##  9    10    10     2
## 10     2     2     6
## # ℹ 17 more rows

Podemos ampliar esto con el parámetro levels. Ahora será con el mínimo, el quantile 33, el cuantil 66 y el máximo.

grid_regular(
  mtry(range  = c(2, 10)),
  min_n(range = c(2, 10)),
  trees(range = c(2, 10)),
  levels = 4  )
## # A tibble: 64 × 3
##     mtry min_n trees
##    <int> <int> <int>
##  1     2     2     2
##  2     4     2     2
##  3     7     2     2
##  4    10     2     2
##  5     2     4     2
##  6     4     4     2
##  7     7     4     2
##  8    10     4     2
##  9     2     7     2
## 10     4     7     2
## # ℹ 54 more rows
grid_max_entropy(
  mtry(range  = c(1, 4)),
  min_n(range = c(10, 30)),
  trees(range = c(1, 1000)),
  size = 10
  )
## # A tibble: 10 × 3
##     mtry min_n trees
##    <int> <int> <int>
##  1     1    25    20
##  2     2    11   610
##  3     4    12     2
##  4     2    20   936
##  5     4    12   733
##  6     4    30   868
##  7     1    20   676
##  8     2    19   307
##  9     2    29   661
## 10     4    26   190

10.2 Ajustando los hiper-parametros: convergen yardstick con dial y tune

WF_hiper_parametros <- 
  # 🔴 Este proceso se hace con la función tune grid  🔴
  tune_grid(
  # Ponemos la receta nuestro modelo, que incuye el preprocesamiento y la fórmula
  object = WF_Receta_modelo_1_Random_forest,
  # 🟡 Ponemos los Folds del conjunto de entrenamiento,  
  # 🟡 Dentro del entrenamiento hará un sub entrenamiento y una sub validación
  resamples = Folds_training,
  #Ponemos un grid para que pruebe con diferentes híper-parámetros, ya con un rango pre-definido
  grid = grid_max_entropy(
    mtry(range  = c(1, 4)),
    min_n(range = c(10, 30)),
    trees(range = c(1, 1000)),
    size = 10),
  # Ponemos las métricas que hemos definido anteriormente
  metrics = Modelo_Metricas, 
  # Con esta opción podemos guardar las predicciones
  control = control_grid( save_pred = T)
)
WF_hiperparametros_coleccion <- 
  WF_hiper_parametros %>% 
  # 🔴 Esta es la funcion que retorna las métricas, importante tenerla en cuenta
  collect_metrics() 

WF_hiperparametros_coleccion %>% 
  head(.,20) %>% 
  knitr::kable()
mtry trees min_n .metric .estimator mean n std_err .config
2 886 10 accuracy binary 0.6610256 5 0.0496335 Preprocessor1_Model01
2 886 10 f_meas binary 0.7673926 5 0.0365573 Preprocessor1_Model01
2 886 10 j_index binary 0.1832661 5 0.1253188 Preprocessor1_Model01
2 886 10 mcc binary 0.1626371 5 0.1086564 Preprocessor1_Model01
2 886 10 precision binary 0.8400397 5 0.0310967 Preprocessor1_Model01
2 886 10 recall binary 0.7082661 5 0.0439273 Preprocessor1_Model01
2 886 10 roc_auc binary 0.7213836 5 0.0469202 Preprocessor1_Model01
2 886 10 sensitivity binary 0.7082661 5 0.0439273 Preprocessor1_Model01
2 886 10 specificity binary 0.4750000 5 0.0918559 Preprocessor1_Model01
4 123 28 accuracy binary 0.6760256 5 0.0501729 Preprocessor1_Model02
4 123 28 f_meas binary 0.7747366 5 0.0350463 Preprocessor1_Model02
4 123 28 j_index binary 0.2768145 5 0.1524681 Preprocessor1_Model02
4 123 28 mcc binary 0.2346557 5 0.1298836 Preprocessor1_Model02
4 123 28 precision binary 0.8684292 5 0.0392609 Preprocessor1_Model02
4 123 28 recall binary 0.7018145 5 0.0393757 Preprocessor1_Model02
4 123 28 roc_auc binary 0.7346900 5 0.0462352 Preprocessor1_Model02
4 123 28 sensitivity binary 0.7018145 5 0.0393757 Preprocessor1_Model02
4 123 28 specificity binary 0.5750000 5 0.1286954 Preprocessor1_Model02
2 81 16 accuracy binary 0.6762821 5 0.0479575 Preprocessor1_Model03
2 81 16 f_meas binary 0.7791279 5 0.0342760 Preprocessor1_Model03
autoplot(WF_hiper_parametros) +
  theme_bw()

WF_hiper_parametros %>% 
  collect_metrics() %>% 
  pivot_longer(cols = c('mtry','trees', 'min_n'), names_to = 'tipo_parametro', values_to = 'valor_parametro') %>% 
  ggplot(aes(valor_parametro,mean, color = mean)) +
  geom_point() +
  facet_grid(.metric~tipo_parametro, scales = 'free' ) +
  scale_color_viridis_b() +
  theme_bw()

WF_hiperparametros_mejor <- WF_hiper_parametros %>%  show_best("roc_auc",1)

WF_hiperparametros_mejor %>% 
  knitr::kable()
mtry trees min_n .metric .estimator mean n std_err .config
2 240 29 roc_auc binary 0.7387727 5 0.0384544 Preprocessor1_Model06

11 La conversión de todo, fit y predict

Ya hemos hecho todos los procesos previos. Todo lo que a un modelo de machine learning le puede hacer falta. Finalmente, una vez que lo tenemos todo hilado, hay que proceder a entrenar el modelo y dejarlo listo.

Modelo_fitted <- WF_Receta_modelo_1_Random_forest %>% 
  # 🔴 Finalizamos el flujo de trabajo con la mejor combinación de Hiper paraemtros encontrada
  finalize_workflow(WF_hiperparametros_mejor) %>%
  # 🔴 finalmente ajustamos con fit
  fit(data = Training_datos)

Con fit ya por fin tenemos nuestro modelo, ya podemos ponernos a hacer lo que hace el aprendizaje automático: predecir. vamos a crear una observación nueva y a ver como la trabaja.

Nuevo_caso <- tibble(
  'Cognicion_01' = 123,
  'Cognicion_02' = 102,
  'Cognicion_03' = 89,
  'Cognicion_04' = 100,
  'Cognicion_05' = 179,
  'Escala_01'    = 12,
  'Escala_02'    = 22,
  'Escala_03'    = 20
)
tibble(
  # 🔴 La función predict de la librería stats tiene sinergias con tidymodels, neecsitra de estos argumentos:
  # - 🟡 El modelo ya ajustado (model_fitted) 
  # - 🟡 Un tibble/dataframe con la nueva prediccion (debe contener las mismas variables)
  # - 🟡 Un tipo de prediccion en los casos de clasificación, si queremos las probabiliadades de cada categoria o la categoría predicha  directamente
  predict(Modelo_fitted, Nuevo_caso, type = "prob" ),
  predict(Modelo_fitted, Nuevo_caso, type = "class")) %>% 
  knitr::kable()
.pred_0 .pred_1 .pred_class
0.9541702 0.0458298 0
tibble(predict(Modelo_fitted, Testing_datos, type = "prob" ), predict(Modelo_fitted, Testing_datos, type = "class")) %>% 
  head(.,10) %>% 
  knitr::kable()
.pred_0 .pred_1 .pred_class
0.7868158 0.2131842 0
0.5619997 0.4380003 0
0.8653084 0.1346916 0
0.4211667 0.5788333 1
0.8627070 0.1372930 0
0.9247541 0.0752459 0
0.0893169 0.9106831 1
0.9033712 0.0966288 0
0.9443273 0.0556727 0
0.7467458 0.2532542 0

11.1 last_fit

Crear las predicciones y todo es un sistema engorroso. Además de que se pierde mucha informaciíon por el camino, ya que todo sigue estando en objetos dispersos. Para ello los desarrolladores crearon last_fit.

La función last_fit permite tener todo unificado en un solo tibble, al tener que ajustarse con un objeto split, esta recibe directamente los datos de entrenamiento y los datos de evaluación (recibe datos de training y testing) y entonces crea sus métricas y sus predicciones.

Con last_fit() contamos con una serie de funciones para extraer todo lo que necesitemos de ellas:

  • extract_fit_engine() : permite extraer el modelo crudo como si no hubiera sido ajustado con tidymodels, de modo que si es un modelo explicativo puede ser tratado con normalidad
  • collect_metrics() : extraer un tibble con las métricas del modelo (puestas con metric set)
  • collect_predictions() : extrae un tibble con las predicciones de los datos de evaluación. (Por defecto el umbral de rendimiento en esta función es 0.5, y estoy buscando alguna forma de cambiar esto, pero se puede manipular el mismo tibble como cualquier otro, de modo que no es un drama.)
  • extract_workflow() : permite extraer todo el flujo del modelo, de incluyendo la receta de modo que podemos hacer predicciones:
Modelo_1_final <-  WF_Receta_modelo_1_Random_forest %>% 
  # 🔴 Finalizamos el flujo de trabajo con la mejor combinación de Hiper paraemtros encontrada
  finalize_workflow(WF_hiperparametros_mejor) %>%
  # 🔴 Ajustamos con last_fit, poniendo el split de datos inicial y las metricas para las predicciones
  last_fit(Split_datos_strat, metrics = Modelo_Metricas)

Modelo_1_final
## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits           id               .metrics .notes   .predictions .workflow 
##   <list>           <chr>            <list>   <list>   <list>       <list>    
## 1 <split [198/87]> train/test split <tibble> <tibble> <tibble>     <workflow>
Modelo_1_final %>%
  # 🔴 dentro del flujo de trabajo están la receta y el modelo ya ajustado
  extract_workflow() %>% 
  # de modo que lo podemos usar para hacer predicciones
  predict(Testing_datos, type = "prob" )
## # A tibble: 87 × 2
##    .pred_0 .pred_1
##      <dbl>   <dbl>
##  1  0.787   0.213 
##  2  0.562   0.438 
##  3  0.865   0.135 
##  4  0.421   0.579 
##  5  0.863   0.137 
##  6  0.925   0.0752
##  7  0.0893  0.911 
##  8  0.903   0.0966
##  9  0.944   0.0557
## 10  0.747   0.253 
## # ℹ 77 more rows

A partir de aquí podemos manipular el objeto fit de múltiples maneras múltiples.

Modelo_1_final %>%
  # 🔴 Con esta funcion podemos recopilar todas las metricas de un modelo
  collect_metrics() %>% 
  select(-.estimator, -.config) %>%  
  mutate(.estimate = round(.estimate,2)) %>% 
  knitr::kable()
.metric .estimate
accuracy 0.71
j_index 0.39
precision 0.89
sensitivity 0.72
specificity 0.67
f_meas 0.80
recall 0.72
mcc 0.33
roc_auc 0.74
options(yardstick.event_first = FALSE)

Modelo_1_final %>%  collect_predictions() %>% 
  # 🔴 Esta funcion permite calcular rapidamente la curva roc de un modelo ya ajustado
  roc_curve(truth = Var_respuesta,  .pred_0 ) %>% 
  ggplot(aes( x = 1 - specificity ,y = sensitivity ) ) +
  geom_path() +
  geom_abline(intercept = 0 , slope = 1, linetype = 3) +
  labs(title = 'Cruva roc del modelo_1_final') +
  theme_bw()

Advertencia: en la función de curva roc, ved que puesto predicción del 0 en vez del 1. Esto es un error interno en la clasificación dentro del del paquete yardstick. Hay que tener cuidado sobre si estas funciones utilizan el primer / segundo / x nivel de su factor como el “evento”. Por defecto, yardstick elige el primer nivel de verdad como “evento” al calcular la curva roc. Podéis encontrar más de este suceso en el github de los desarrolladores: https://github.com/tidymodels/yardstick/issues/94

12 Rendimiento del umbral (threshold performance en modelos de clasificación)

El concepto de rendimiento del umbral (threshold performance) es crucial a la hora de evaluar y ajustar modelos de clasificación binaria.

El umbral es el valor de probabilidad por encima del cual un resultado predicho se considera de clase positiva. Por defecto, muchos modelos utilizan un umbral de 0,5. Sin embargo, en casi todos los casos, el umbral óptimo no estará otro punto. Este “óptimo lo tendremos que definir que en función de los objetivos y requisitos específicos del problema.

Ajustar el umbral permite elegir entre precisión y la recall. Bajar el umbral aumenta la sensibilidad (recall) pero disminuye la precisión, y viceversa. Básicamente, el umbral óptimo depende de la importancia relativa de los falsos positivos y los falsos negativos en una aplicación concreta.

En tidymodels está la librería probably con la función threshold_perf(), la cual nos permite, dadas las predicciones del modelo ver como se calcula este umbral.

pacman::p_load(probably)


Modelo_1_final_Threshold_performance <- 
  # Llamamos al modelo ya finallizado con "last_fit"
  Modelo_1_final %>% 
  # Coleccionamos sus predicciones
  collect_predictions() %>% 
  # 🔴 Ejecutamos el redimeinto de umbral con la siguiente función
  threshold_perf(
    # Le decimos que variable es la verdadera respuesta
    truth = Var_respuesta,  
    # le decimos que variable es la predicción (aviso: debido a un bugg en la version )
    .pred_1, 
    # ajustamos un rango en el que se evaluará el umbral en cada punto
    thresholds = seq(0.5,1,0.01) )

# displayment
Modelo_1_final_Threshold_performance %>% 
  head(., 10) %>% 
  mutate(.estimate = round(.estimate, 2)) %>% 
  knitr::kable()
.threshold .metric .estimator .estimate
0.50 sensitivity binary 0.28
0.51 sensitivity binary 0.28
0.52 sensitivity binary 0.28
0.53 sensitivity binary 0.28
0.54 sensitivity binary 0.28
0.55 sensitivity binary 0.28
0.56 sensitivity binary 0.26
0.57 sensitivity binary 0.25
0.58 sensitivity binary 0.23
0.59 sensitivity binary 0.23
Modelo_1_final_Threshold_performance %>% 
  ggplot(aes(x = .threshold, y = .estimate, color = .metric)) +
  geom_line(size = 1.1) +
  scale_color_manual(values = c('purple','blue','pink')) +
  geom_vline(xintercept = 
    Modelo_1_final_Threshold_performance %>%  
      filter(.metric == 'j_index') %>%
      arrange(desc(.estimate)) %>% 
      slice(1) %>% 
      pull(.threshold), linetype = 'dashed', size = 1.5)

13 El cielo es el límite, workflow_set()

Bien, ya hemos aprendido a crear todo un flujo entero de machine learning con tidymodels, desde la base hasta su final. Muy útil todo. pero calibrar modelos de 1 en 1 es un auténtico petardazo. No me acaba de ser útil este taller. ¿Cómo me va a ser útil si tengo que probar unas 20 fórmulas con datos?

Lo entiendo perfectamente, a mi también me pasó. Los médicos me preguntaron que era mejor para tener en cuenta, si las escalas, las cogniciones, la bioquímica o la genética a la hora de diagnosticar al paciente.

Comprobémoslo…

13.1 Todas las recetas

Receta_modelo_Escalas_formula <- Training_datos %>% 
  dplyr::select(dplyr::matches('Esc')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Escalas <- recipe(  
    formula = Receta_modelo_Escalas_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_bag(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()



Receta_modelo_Cognicion_formula <- Training_datos %>% 
  dplyr::select(
   dplyr::matches('Cog')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Cognicion <- recipe(  
    formula = Receta_modelo_Cognicion_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>%
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_mean(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

bake(Receta_modelo_Cognicion, new_data = NULL)
## # A tibble: 316 × 4
##    Cognicion_01 Cognicion_03 Cognicion_04 Var_respuesta
##           <dbl>        <dbl>        <dbl> <fct>        
##  1       -0.643    -2.70e-15     2.51e-15 0            
##  2        0.760     3.13e- 1     1.13e+ 0 0            
##  3       -2.91     -2.70e-15     2.51e-15 0            
##  4       -1.11     -2.70e-15     2.51e-15 0            
##  5        0.797     1.64e+ 0     2.59e- 1 0            
##  6        0.673     5.70e- 1     3.93e- 1 0            
##  7        0.203    -2.70e-15     2.51e-15 0            
##  8        0.211     2.03e- 1     8.86e- 1 0            
##  9       -0.545    -2.70e-15     2.51e-15 0            
## 10        0.873     1.50e+ 0     1.09e+ 0 0            
## # ℹ 306 more rows
Receta_modelo_Bioquimicas_formula <- Training_datos %>% 
  dplyr::select(
   dplyr::matches('Bio')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Bioquimica <- recipe(  
    formula = Receta_modelo_Bioquimicas_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_mean(all_numeric()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() )  %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

Receta_modelo_Popurri <- Training_datos %>%
  dplyr::select(
    Cognicion_01, Cognicion_02, Escala_01,Escala_02, Bioquimica_01,Bioquimica_02, Genetica_01,Genetica_02) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Popurri <- recipe(  
    formula = Receta_modelo_Popurri,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 5% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_bag(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

13.2 workflow_set() para unificarlas todas

Ahora, esto son muchas recetas y mucho texto. Será difícil combinarlo todo, ¿no?…¡Qué va!. workflow_set() ya está pensado para esto. Funciona como el crossing() de la librería purrr o como el expand.grid() de la librería base, crear workflows con todas las combinaciones.

Receta_list <- list(
  'Modelo_Escalas'    = Receta_modelo_Escalas,
  'Modelo_Cognicion'  = Receta_modelo_Cognicion,
  'Modelo_Bioquimica' = Receta_modelo_Bioquimica,
  'Modelo_Popurri'    = Receta_modelo_Popurri
  )
  
Model_list <- list(
  RandomForest     = Model_RandomForest
  # ,XGBoost          = Model_XGBoost
  # ,Bagged_Trees     = Model_Bagged_Tree
  )

Wokflows_set <- workflow_set(
  preproc = Receta_list, 
  models = Model_list, 
  # 🔴 La opción Cross es para que haga todos los cruces de modelos con recetas
  cross = T)


Wokflows_set
## # A workflow set/tibble: 4 × 4
##   wflow_id                       info             option    result    
##   <chr>                          <list>           <list>    <list>    
## 1 Modelo_Escalas_RandomForest    <tibble [1 × 4]> <opts[0]> <list [0]>
## 2 Modelo_Cognicion_RandomForest  <tibble [1 × 4]> <opts[0]> <list [0]>
## 3 Modelo_Bioquimica_RandomForest <tibble [1 × 4]> <opts[0]> <list [0]>
## 4 Modelo_Popurri_RandomForest    <tibble [1 × 4]> <opts[0]> <list [0]>

Una vez hecho el workflow_set, se pone a trabajar con workflow_map.

Wokflows_set_map <- 
  Wokflows_set %>% 
  # 🔴 workflow maps es una funcion que permite ajustar múltiples flujos de trabajo
    workflow_map(
    resamples = Folds_training, 
    fn = "tune_grid",
    grid = grid_max_entropy(
      mtry(range  = c(1, 4)),
      min_n(range = c(10, 30)),
      trees(range = c(1, 1000)),
      size = 10),
    # verbose = TRUE, 
    metrics = Modelo_Metricas, 
    control = control_grid( save_pred = T),
    # 🔴 para garantizare replicabildiad, podemos fijar la semilla en el proceso
    seed = 2465)

Wokflows_set_map
## # A workflow set/tibble: 4 × 4
##   wflow_id                       info             option    result   
##   <chr>                          <list>           <list>    <list>   
## 1 Modelo_Escalas_RandomForest    <tibble [1 × 4]> <opts[4]> <tune[+]>
## 2 Modelo_Cognicion_RandomForest  <tibble [1 × 4]> <opts[4]> <tune[+]>
## 3 Modelo_Bioquimica_RandomForest <tibble [1 × 4]> <opts[4]> <tune[+]>
## 4 Modelo_Popurri_RandomForest    <tibble [1 × 4]> <opts[4]> <tune[+]>

AVISO, workflow_map puede paralelizarse, pero aún estoy descubriendo como integrarlo de forma efectiva.

Finalmente, podemos calibrar cada modelo individualmente y tener todo finiquitado. Ese proceso sería más eficiente con una función, desde luego. Pero se me acaba el tiempo para dedicar al taller :(

# Modelo Escalas ----
Modelo_Escalas_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  # 🔴 Extraemos el flujo correspondiente al modelo que queremos 
  extract_workflow_set_result('Modelo_Escalas_RandomForest') %>%
  # 🔴 podemos ver una tabla con los mejores resultados, acorde a una metrica
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry,trees,min_n,.config )

Modelo_Escalas_RandomForest_final <- Wokflows_set_map %>% 
  # 🔴 Extraemos el flujo correspondiente al modelo que queremos 
  extract_workflow('Modelo_Escalas_RandomForest') %>% 
  finalize_workflow(Modelo_Escalas_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

# Modelo_Cognicion ----

Modelo_Cognicion_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  extract_workflow_set_result('Modelo_Cognicion_RandomForest') %>%
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry, trees, min_n, .config )

Modelo_Cognicion_RandomForest_final <- Wokflows_set_map %>% 
  extract_workflow('Modelo_Cognicion_RandomForest') %>% 
  finalize_workflow(Modelo_Cognicion_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

# Modelo_Bioquimico ----

Modelo_Bioquimica_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  extract_workflow_set_result('Modelo_Bioquimica_RandomForest') %>%
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry, trees, min_n, .config )

Modelo_Bioquimica_RandomForest_final <- Wokflows_set_map %>% 
  extract_workflow('Modelo_Bioquimica_RandomForest') %>% 
  finalize_workflow(Modelo_Bioquimica_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

# Modelo Popurri ----
Modelo_Popurri_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  extract_workflow_set_result('Modelo_Popurri_RandomForest') %>%
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry, trees, min_n, .config )

Modelo_Popurri_RandomForest_final <- Wokflows_set_map %>% 
  extract_workflow('Modelo_Popurri_RandomForest') %>% 
  finalize_workflow(Modelo_Popurri_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

Una vez acabado todo, podemos ver todo lo necesario usando unos pocos verbos como collect() y las lógicas del tidyverse.

13.3 Métricas de todos los modelos

map2(
  list(
    Modelo_Escalas_RandomForest_final,
    Modelo_Cognicion_RandomForest_final,
    Modelo_Bioquimica_RandomForest_final,
    Modelo_Popurri_RandomForest_final),
  list('Modelo_Escala','Modelo_Cognicion','Modelo_Bioquimica','Modelo_Popurri'),
  ~ ..1 %>% 
    collect_metrics() %>% 
    mutate('Modelo'= ..2) %>% 
    select(.metric,.estimate,Modelo )) %>% 
  bind_rows() %>% 
  pivot_wider(names_from = Modelo, values_from = .estimate) %>% 
  knitr::kable()
.metric Modelo_Escala Modelo_Cognicion Modelo_Bioquimica Modelo_Popurri
accuracy 0.6436782 0.6091954 0.7471264 0.6666667
j_index 0.2222222 -0.0265700 0.0241546 0.2101449
precision 0.8518519 0.7868852 0.7974684 0.8448276
sensitivity 0.6666667 0.6956522 0.9130435 0.7101449
specificity 0.5555556 0.2777778 0.1111111 0.5000000
f_meas 0.7479675 0.7384615 0.8513514 0.7716535
recall 0.6666667 0.6956522 0.9130435 0.7101449
mcc 0.1855216 -0.0235126 0.0338612 0.1805788
roc_auc 0.7504026 0.5024155 0.5048309 0.6231884

13.4 Cruvas ROC todos los modelos

ROC_CURVE_training_todos_los_posibles_modelos <- Wokflows_set_map %>%  
  collect_predictions() %>%
  group_split(wflow_id) %>% 
  set_names(list('Escala','Cognicion','Bioquimica','Popurri')) %>% 
  map(~ .x %>% roc_curve(truth = Var_respuesta , .pred_0)) %>% 
  map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
               ~ .x %>%  mutate(model = .y)) %>% 
         bind_rows() 

Plot_ROC_CURVE_training_todos_los_posibles_modelos <- ROC_CURVE_training_todos_los_posibles_modelos %>% 
  ggplot(aes( x = 1 - specificity,y = sensitivity, color = model ) ) +
  geom_path() +
  geom_abline(intercept = 0, slope = 1, linetype = 3) +
  theme_bw()
  


ROC_CURVE_training_modelos_finales <- list(
  list('Esc' ,Modelo_Escalas_RandomForest_hyperparametros$.config),
  list('Cog' ,Modelo_Cognicion_RandomForest_hyperparametros$.config),
  list('Bio' ,Modelo_Bioquimica_RandomForest_hyperparametros$.config),
  list('Pop' ,Modelo_Popurri_RandomForest_hyperparametros$.config)) %>% 
  map(~ Wokflows_set_map %>% 
        collect_predictions() %>% 
        filter(
          str_detect(wflow_id, .x[[1]]) & 
            .config== .x[[2]]
        )) %>% 
  
  map(~ .x %>%
        mutate(Model= str_extract(wflow_id, '^.{2}') ) %>% 
        select(Model,Var_respuesta, .pred_0, .pred_1, .pred_class)
  ) %>% 
  set_names(list('Escala','Cognicion','Bioquimica','Popurri') ) %>%
  map(~ .x %>% roc_curve(truth = Var_respuesta, .pred_0)) %>% 
  map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
       ~ .x %>% mutate(model= .y)) %>% 
  bind_rows() 
  
Plot_ROC_CURVE_training_modelos_finales <- ROC_CURVE_training_modelos_finales %>% 
  ggplot(aes( x = 1 - specificity ,y = sensitivity, color = model ) ) +
  geom_path() +
  geom_abline(intercept = 0, slope = 1, linetype = 3) +
  theme_bw()


ROC_CURVE_testing_modelos_finales <- list(
  Modelo_Escalas_RandomForest_final,
  Modelo_Cognicion_RandomForest_final,
  Modelo_Bioquimica_RandomForest_final,
  Modelo_Popurri_RandomForest_final
  ) %>% 
  map(~.x %>% 
        collect_predictions() %>% 
        roc_curve(Var_respuesta, .pred_0)) %>% 
  map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
       ~ .x %>%  mutate(model = .y)) %>% 
  bind_rows()

Plot_ROC_CURVE_testing_modelos_finales <- ROC_CURVE_testing_modelos_finales %>% 
  ggplot(aes(x = 1 - specificity ,y = sensitivity, color = model ) ) +
  geom_path() +
  geom_abline(intercept = 0, slope = 1, linetype = 3) +
  theme_bw()
ggpubr::ggarrange( Plot_ROC_CURVE_training_todos_los_posibles_modelos,
                   Plot_ROC_CURVE_training_modelos_finales, 
                   Plot_ROC_CURVE_testing_modelos_finales ,
                   nrow = 1 , 
                   common.legend = T,
                   legend = "bottom",
                   labels = c('Todos_training','Training','Test')) 

14 Haciendo el machine learning entendible: SHAP Values

Los valores de Shapley (shap values), que deben su nombre al premio Nobel Lloyd Shapley, son un concepto de la teoría de juegos cooperativos que ha acabado encontrado aplicaciones en el aprendizaje automático.

En teoría de juegos, los valores de Shapley son la distribución equitativa de la contribución de cada jugador en un juego cooperativo entre todos los jugadores. En el contexto del aprendizaje automático, las características o variables se consideran jugadores, y el juego consiste en cuánto contribuye cada característica a la predicción de un modelo.

Es decir, los valores de Shapley son sencillamente las marginales de cada variable dentro del modelo. Compilarlos es computacionalmente extenso, ya que hay que hacer cada combinación de valores para una fila y evaluar de cuanto es su marginal. El valor final llamado Shap value es promedio de todas las aportaciones de una variable marginal.

Aparte de su extenso coste computacional, todo shap value asume que está utilizando un modelo de caja negra. de manera que si el modelo no está bien especificado en tema de interacciones, etc. es posible que estén dando indicaciones erróneas.

Actualmente se pueden compilar en tidymodels utilizando la librería agua, al ser una extensión de tidymodels del entorno de h2o. El problema es que esto hace que el motor de los modelos deba ser a la fuerza h2o.

Una propuesta sería buscar un modelo bien calibrado con tidymodels y pasarlo por h2o para hacer un diagnostico más eficaz.

h2o_start()
h2o::h2o.no_progress()

modelo_h20 <- rand_forest() %>% 
  set_engine("h2o", max_runtime_secs = 20) %>% 
  set_mode('classification') %>% 
  fit(Var_respuesta ~ ., data = Receta_modelo_Bioquimica %>% bake(new_data = NULL) )

entrada_h20 <- h2o::as.h2o(
  Receta_modelo_Bioquimica %>% bake(new_data = NULL) %>%
    dplyr::rename(.outcome = Var_respuesta))

modelo_h20 %>%  
  extract_fit_engine() %>% 
  h2o::h2o.shap_summary_plot(entrada_h20)

# modelo_h20 %>%  
#   extract_fit_engine() %>% 
#   h2o::h2o.explain(entrada_h20)

h2o_end()

Para cada predicción, los valores de Shapley representan la contribución de cada característica a la diferencia entre la predicción del modelo y la predicción media. En otras palabras, indican cuánto añade o resta cada característica a la predicción media.

Si una predicción media es de 0.5, y el shap value es de 0.1, indicará que esa variable sola podría llegar a tener un impacto tal que haga que la predicción llegue a 0.6.